ltbl <- NULL
ltbl[1:3] <- "Reproduction driven expansion"
ltbl[4:6] <- "Death driven expansion"

# Set ylim on divergence
ylimit <- 0.2
xlimit <- 40
palette(c("red", "green", "blue"))
par(mfrow=c(2,5))
for(i in 1:1){
	divergence <- read.csv(paste("sweep",i,"/sweep",i,".divergence.log",sep=""),header=FALSE)
	biopsy.est.diverg <- read.csv(paste("sweep",i,"/sweep",i,".biopsy.divergence.log",sep=""),header=FALSE)

	for(biop in 1:max(biopsy.est.diverg[,2])){

	plot(0,0,t="n",xlab="Number of Cells",ylab="Divergence",ylim=c(0,ylimit),xlim=c(0,xlimit))
	title(main=ltbl[i],line=1)
	for(repl in 1:max(biopsy.est.diverg[,1])){
	   lin <- biopsy.est.diverg[biopsy.est.diverg[,2]==biop&biopsy.est.diverg[,1]==repl,]
	   points(lin[,3],lin[,4],col="black",xlab="Number of Cells",ylab="Divergence",ylim=c(0,ylimit),xlim=c(0,xlimit))
	   abline(h=divergence[divergence[,1]==7300,2], col="red")
	}
	}
}


# Set ylim on divergence
ylimit <- 0.15
xlimit <- 5
palette(c("red", "green", "blue"))
par(mfrow=c(2,5))
for(i in 1:1){
	divergence <- read.csv(paste("sweep",i,"/sweep",i,".divergence.log",sep=""),header=FALSE)
	biopsy.est.diverg <- read.csv(paste("sweep",i,"/sweep",i,".biopsy.divergence.log",sep=""),header=FALSE)

	for(biop in 1:max(biopsy.est.diverg[,2])){

	plot(0,0,t="n",xlab="Number of Cells",ylab="Divergence",ylim=c(0,ylimit),xlim=c(0,xlimit))
	title(main=ltbl[i],line=1)
	medianslin <- NULL
	lowerlin <- NULL
	upperlin <- NULL

	for(cell in 1:max(biopsy.est.diverg[,3])){
	   lin <- biopsy.est.diverg[biopsy.est.diverg[,2]==biop&biopsy.est.diverg[,3]==cell,]
	   medianslin[cell] <- median(lin[,4])
	   lowerlin[cell] <- as.numeric(quantile(lin[,4])[1])
	   upperlin[cell] <- as.numeric(quantile(lin[,4])[4])
	}
	   linx <- c(1:max(biopsy.est.diverg[,3]))

	   lines(linx,medianslin,lty=1,col="black")
	   lines(linx,lowerlin,lty=2,col="black")
	   lines(linx,upperlin,lty=2,col="black")
	   abline(h=divergence[divergence[,1]==7300,2], col="red")

	}
}


ltbl <- NULL
ltbl[1:3] <- "Reproduction driven expansion"
ltbl[4:6] <- "Death driven expansion"
# Set ylim on divergence
ylimit <- 0.70
xlimit <- 5
palette(c("red", "green", "blue"))
par(mfrow=c(2,5))
for(i in 4:4){
	divergence <- read.csv(paste("sweep",i,"/sweep",i,".divergence.log",sep=""),header=FALSE)
	biopsy.est.diverg <- read.csv(paste("sweep",i,"/sweep",i,".cutout.biopsy",sep=""),header=FALSE)

	for(biop in 1:5){

	plot(0,0,t="n",xlab="Number of Cells",ylab="Divergence",ylim=c(0,ylimit),xlim=c(0,xlimit))
	title(main=ltbl[i],line=1)
	medianslin <- NULL
	lowerlin <- NULL
	upperlin <- NULL

	for(cell in 1:max(biopsy.est.diverg[,3])){
	   lin <- biopsy.est.diverg[biopsy.est.diverg[,2]==biop&biopsy.est.diverg[,3]==cell,]
	   medianslin[cell] <- median(lin[,4])
	   lowerlin[cell] <- as.numeric(quantile(lin[,4])[1])
	   upperlin[cell] <- as.numeric(quantile(lin[,4])[4])
	}
	   linx <- c(1:max(biopsy.est.diverg[,3]))

	   lines(linx,medianslin,lty=1,col="black")
	   lines(linx,lowerlin,lty=2,col="black")
	   lines(linx,upperlin,lty=2,col="black")
	   abline(h=divergence[divergence[,1]==7300,2], col="red")

	}
}









# Plot clone size distribution over time
par(mfrow=c(1,1))
for(i in 2:2){
	sizedistr <- read.csv(paste("sweep",i,"/sweep",i,".unique.clones.size.distrib.log",sep=""))
	boxplot(split(sizedistr[,2], sizedistr[,1]),ylim=c(0,200))
	clonenums <- as.numeric(sapply(split(sizedistr[,2], sizedistr[,1]),length))
}








# Set ylim on divergence
ylimit <- 0.04
# s=0.5, selmu = 10^-6, 10^-7, 10^-8 # s=1, selmu = 10^-6, 10^-7, 10^-8
par(mfrow=c(2,3))
for(i in 7:8){
	divergence <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".divergence.log",sep=""))
	plot(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="black")
	fixations <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="black",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="red")
	fixations <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="red",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="blue")
	fixations <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="blue",pch=fixev[,5]+14,cex=1.5)
}

# Repro driven
# Set ylim on divergence
ylimit <- 0.04
# s=0.5, selmu = 10^-6, 10^-7, 10^-8 # s=1, selmu = 10^-6, 10^-7, 10^-8
par(mfrow=c(2,3))
for(i in 9:12){
	divergence <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".divergence.log",sep=""))
	plot(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="black")
	fixations <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".fixations.log",sep=""))
	divergence <- divergence[1:nrow(fixations),]
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="black",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="red")
	fixations <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".fixations.log",sep=""))
	divergence <- divergence[1:nrow(fixations),]
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="red",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="blue")
	fixations <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".fixations.log",sep=""))
	divergence <- divergence[1:nrow(fixations),]
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="blue",pch=fixev[,5]+14,cex=1.5)
}

# Repro driven
# Set ylim on divergence
ylimit <- 0.04
# s=0.5, selmu = 10^-6, 10^-7, 10^-8 # s=1, selmu = 10^-6, 10^-7, 10^-8
#par(mfrow=c(2,3))
for(i in 15:15){
	divergence <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".divergence.log",sep=""))
	plot(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="black")
	fixations <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="black",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="red")
	fixations <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="red",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="blue")
	fixations <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="blue",pch=fixev[,5]+14,cex=1.5)
}









for(i in 3:5){
	divergence <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".divergence.log",sep=""))
	plot(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="black")
	fixations <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".fixations.log",sep=""))
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],divergence[,2],col=fixev[,2])

	divergence <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="red")
	fixations <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".fixations.log",sep=""))
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],divergence[,2],col=fixev[,2])

	divergence <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="blue")
	fixations <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".fixations.log",sep=""))
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],divergence[,2],col=fixev[,2])
}




d <- read.csv("sweep21/sweep21.frequencies",header=FALSE)

dat <- d[d[,3]>0.05&d[,1]>0,]

rrr <- unique(dat[,1])
rrrc <- gray(rrr/max(rrr))
clonecol <- cbind(rrr,rrrc)

plot(dat[,2],dat[,3],col="black",ylab="Frequency",xlab="Time (days)",pch=19)
for(i in 1:length(rrr)){
lines(dat[dat[,1]==rrr[i],2],dat[dat[,1]==rrr[i],3],col="black")
}

# Color by clone selective advantage

# Load clones selective advantage, color clones by advantage
clones.fitness <- read.csv("sweep21/sweep21.clones",header=FALSE)
d <- read.csv("sweep21/sweep21.frequencies",header=FALSE)

# Remove duplicates
clones.fitness <- unique(clones.fitness)

cfit <- cbind(rep(0,nrow(clones.fitness)),rep(0,nrow(clones.fitness)))
cfit[clones.fitness[,2]==1,1] <- "black"
cfit[clones.fitness[,2]==2,1] <- "orange"
cfit[clones.fitness[,2]==4,1] <- "red"
cfit[clones.fitness[,2]==8,1] <- "brown"
cfit[clones.fitness[,3]==1,2] <- "black"
cfit[clones.fitness[,3]==2,2] <- "orange"
cfit[clones.fitness[,3]==4,2] <- "red"
cfit[clones.fitness[,3]==8,2] <- "brown"

# get clones that have reached more than 5% frequency over 20 years
# clone_id, time, frequency
dat <- d[d[,3]>0.05&d[,1]>0,]

intersect <- function(x, y) y[match(x, y, nomatch = 0)]

# get the frequency of clones in dat, even if it dropped to 0 at some time points
clones.freq <- d[d[,1] %in% dat[,1],]

# Plot
par(mfrow=c(1,2))
plot(clones.freq[,2],clones.freq[,3],col="black",ylab="Frequency",xlab="Time (days)",main="Survival Advantage",pch=NA)
rrr <- unique(clones.freq[,1])
for(i in 1:length(rrr)){
lines(clones.freq[clones.freq[,1]==rrr[i],2],clones.freq[clones.freq[,1]==rrr[i],3],col=cfit[clones.fitness[,1]==rrr[i],1])
}
plot(clones.freq[,2],clones.freq[,3],col="black",ylab="Frequency",xlab="Time (days)",main="Reproductive Advantage",pch=NA)
rrr <- unique(clones.freq[,1])
for(i in 1:length(rrr)){
lines(clones.freq[clones.freq[,1]==rrr[i],2],clones.freq[clones.freq[,1]==rrr[i],3],col=cfit[clones.fitness[,1]==rrr[i],2])
}



# Load data
clonesn21 <- read.csv("sweep21/sweep21.clones.number.log",header=FALSE)

par(mfrow=c(1,1))
plot(clonesn21[,1],clonesn21[,2],col="black",ylab="Absolute number of clones",xlab="Time (days)",main="10e-06",type="l")




# Plot number of clones over time
# Load data
clonesn6 <- read.csv("sweep3/sweep3.clones.number.log",header=FALSE)
clonesn7 <- read.csv("sweep5/sweep5.clones.number.log",header=FALSE)
clonesn8 <- read.csv("sweep6/sweep6.clones.number.log",header=FALSE)

par(mfrow=c(1,3))
plot(clonesn6[,1],clonesn6[,2],col="black",ylab="Absolute number of clones",xlab="Time (days)",main="10e-06",type="l")
plot(clonesn7[,1],clonesn7[,2],col="black",ylab="Absolute number of clones",xlab="Time (days)",main="10e-07",type="l")
plot(clonesn8[,1],clonesn8[,2],col="black",ylab="Absolute number of clones",xlab="Time (days)",main="10e-08",type="l")

grid6 <- read.csv("sweep6/sweep6.log",header=FALSE)
ggrid6 <- grid6[,2:3]
cnum <- NULL
ctime <- NULL
for(i in 1:40){
	cnum[i] <- nrow(unique(ggrid6[ggrid6[,2]==i*182,1:2]))
	ctime[i] <- i*182
}
par(mfrow=c(1,1))
plot(ctime,cnum,col="black",ylab="Absolute number of clones in grid",xlab="Time (days)",main="10e-08",type="l")







# Average cell lifespan
distrib <- NULL
for(i in 1:1000){

	lifespan <- 0	
	div <- rexp(1,1/4)
	death <- rexp(1,1/40)
	while(div<death){
		lifespan <- lifespan + 1
		div <- rexp(1,1/4)
		death <- rexp(1,1/40)
	}	
	distrib[i] <- lifespan
}
hist(distrib,freq=FALSE)






